///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file pfaffian.cpp
///
/// \brief   contains implementations of functions related to pfaffian calculations for the
///          generalized pairing wavefunction
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include <src/pcps.h>
#include <src/blas_lapack.h>
#include <src/matrix_algebra.h>
#include <src/pfaffian.h>

#include <pfapack.h> // interface to Wimmer's pfaffian computation library

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Populates row_indices, which holds each site's row index in the occupied pairing
///          matrix (unoccupied sites are marked with -1).  Also populates occ_sites, which
///          holds the site corresponding to each row of the pairing matrix.
///
/// \param[in]     np           number of particles
/// \param[in]     ns           number of lattice sites
/// \param[in]     config       lattice configuration, size ns
/// \param[out]    occ_sites    on exit, the list of occupied sites, size np
/// \param[out]    row_indices  on exit, the row indices for occupied sites, size ns
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void pcps::pfn::init_occ_sites(const int np,
                               const int ns,
                               const int * const config,
                               int * const occ_sites,
                               int * const row_indices)
{
  int count = 0;
  for (int i = 0; i < ns; i++) {
    if (config[i] == 1) {
      occ_sites[count] = i;
      row_indices[i] = count;
      count++;
    } else {
      row_indices[i] = -1;
    }
  }
  assert(count == np);
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Builds the occupied pairing matrix for the given set of occupied sites.
///
/// \param[in]     np           number of particles
/// \param[in]     ns           number of lattice sites
/// \param[in]     pairing_mat  the full skew-symmetric pairing matrix, size ns*ns
/// \param[in]     occ_sites    the list of occupied sites, size np
/// \param[out]    mat          on exit, the occupied pairing matrix, size np*np
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::pfn::build_occ_par_mat(const int np,
                                                          const int ns,
                                                          const S * const pairing_mat,
                                                          const int * const occ_sites,
                                                          S * const mat)
{
  for (int i = 0; i < np; i++)
  for (int j = 0; j < np; j++)
    mat[i*np+j] = pairing_mat[ occ_sites[i] * ns + occ_sites[j] ];
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Exchanges the rth row and column of the inverse of the occupied pairing matrix
///          with the first row and column.
///
/// \param[in]     np           number of particles
/// \param[in]     ns           number of lattice sites
/// \param[in]     r            index of the row and column to exchange with first row and column
/// \param[in,out] inv_mat      Size np*np.
///                             On entry, the inverse of the occupied pairing matrix.
///                             On exit, the 1st and rth rows and columns have been exchanged.
/// \param[in,out] occ_sites    Size np.
///                             On entry, the list of occupied sites.
///                             On exit, the updated list that reflects the row/column exchanges.
/// \param[in,out] row_indices  Size ns.
///                             On entry, the occupied sites' row indices.
///                             On exit, the updated list that reflects the row/column exchanges.
/// \param[out]    wv1          Work space, size np.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::pfn::make_first_row_col(const int np,
                                                           const int ns,
                                                           const int r,
                                                           S * const inv_mat,
                                                           int * const occ_sites,
                                                           int * const row_indices,
                                                           S * const wv1)
{

  // we only need to take action if the requested row is not the first row
  if (r != 0) {

    // swap rows
    pcps::xcopy(np, inv_mat + r*np, 1,            wv1, 1);
    pcps::xcopy(np, inv_mat + 0*np, 1, inv_mat + r*np, 1);
    pcps::xcopy(np,            wv1, 1, inv_mat + 0*np, 1);

    // swap columns
    pcps::xcopy(np, inv_mat + r, np,         wv1,  1);
    pcps::xcopy(np, inv_mat + 0, np, inv_mat + r, np);
    pcps::xcopy(np,         wv1,  1, inv_mat + 0, np);

    // update indexing
    std::swap(row_indices[occ_sites[0]], row_indices[occ_sites[r]]);
    std::swap(occ_sites[0], occ_sites[r]);

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   For two lattice configurations that differ at the specified sites, this function
///          determines which sites are newly occupied and which are newly empty.
///
/// \param[in]     np             number of particles
/// \param[in]     ns             number of lattice sites
/// \param[in]     nd             the number of sites that differ between the configurations
/// \param[in]     diff_sites     the indices of the differing sites ( 0 <= index < ns ), size nd
/// \param[in]     new_config     new lattice configuration, size ns
/// \param[in]     old_config     old lattice configuration, size ns
/// \param[out]    new_occ_sites  On exit, the newly occupied lattice sites, size 2.
/// \param[out]    old_occ_sites  On exit, the newly empty lattice sites, size 2.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void pcps::pfn::get_diff_occs(const int np,
                              const int ns,
                              const int nd,
                              const int * const diff_sites,
                              const int * const new_config,
                              const int * const old_config,
                              int * const new_occ_sites,
                              int * const old_occ_sites)
{
  assert(nd == 2);
  new_occ_sites[0] = ( old_config[diff_sites[0]] == 0 ? diff_sites[0] : diff_sites[1] );
  old_occ_sites[0] = ( new_config[diff_sites[0]] == 0 ? diff_sites[0] : diff_sites[1] );
  assert( new_config[new_occ_sites[0]] == 1 && old_config[new_occ_sites[0]] == 0 );
  assert( new_config[old_occ_sites[0]] == 0 && old_config[old_occ_sites[0]] == 1 );
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Initializes the inverse of the occupied pairing matrix for the given configuration.
///          Also populates the occupied site list and the list of occupied sites' row indices.
///
/// \param[in]     np             number of particles
/// \param[in]     ns             number of lattice sites
/// \param[in]     config         lattice configuration, size ns   
/// \param[in]     pairing_mat    the full skew-symmetric pairing matrix, size ns*ns
/// \param[out]    occ_sites      on exit, the list of occupied sites, size np   
/// \param[out]    row_indices    on exit, the row indices for occupied sites, size ns   
/// \param[out]    inv_mat        on exit, the inverse of the occupied pairing matrix, size np*np
/// \param[out]    work           work space used in computing inv_mat, size np*np
/// \param[out]    ipiv           work space for pivoting during the LU decomposition, size np   
/// \param[out]    inv_perm       work space for inverting the LU pivoting, size np   
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::pfn::init_matrices(const int np,
                                                      const int ns,
                                                      const int * const config,
                                                      const S * const pairing_mat,
                                                      int * const occ_sites,
                                                      int * const row_indices,
                                                      S * const inv_mat,
                                                      S * const work,
                                                      int * const ipiv,
                                                      int * const inv_perm)
{

  // Populate row_indices and occ_sites
  pcps::pfn::init_occ_sites(np, ns, config, occ_sites, row_indices);

  // build the occupied pairing matrix
  pcps::pfn::build_occ_par_mat(np, ns, pairing_mat, occ_sites, inv_mat);

  // compute the inverse matrix and determinant using an LU decomposition
  S unused;
  pcps::matrix_inverse_lu(np, unused, inv_mat, work, ipiv, inv_perm);

}

// explicitly instantiate the function
template void pcps::pfn::init_matrices(const int np,
                                       const int ns,
                                       const int * const config,
                                       const double * const pairing_mat,
                                       int * const occ_sites,
                                       int * const row_indices,
                                       double * const inv_mat,
                                       double * const work,
                                       int * const ipiv,
                                       int * const inv_perm);
template void pcps::pfn::init_matrices(const int np,
                                       const int ns,
                                       const int * const config,
                                       const std::complex<double> * const pairing_mat,
                                       int * const occ_sites,
                                       int * const row_indices,
                                       std::complex<double> * const inv_mat,
                                       std::complex<double> * const work,
                                       int * const ipiv,
                                       int * const inv_perm);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Computes the pfaffian of a real matrix A.
///
/// \param[in]     n        Dimension of the matrix.
/// \param[in,out] A        On entry, the n*n matrix whose pfaffian is to be computed, which is
///                         assumed to be skew-symmetric.  On exit, the matrix is overwritten.
/// \param[out]    iwork    Integer workspace, size n.
///
/// \return The pfaffian of A in mantissa/exponent format.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
pcps::MantExp<double> pcps::pfn::compute_pfaffian(const int n,
                                                  double * const A,
                                                  int * const iwork)
{

  // odd dimension matrices have pfaffian equal to zero
  if (n % 2 != 0)
    return pcps::MantExp<double>(0.0);

  // compute the pfaffian
  int info;
  double pfaffian[2];
  double work[2];
  int lwork = 2;
  dskpf10_("U", "P", &n, A, &n, pfaffian, iwork, work, &lwork, &info);

  // if the pfaffian fails to compute, return zero so that initial configurations using the pfaffian value will be rejected
  if (info != 0) {
    return pcps::MantExp<double>(0.0);
    //throw pcps::Exception("dskpf10_ failed in pcps::pfn::compute_pfaffian with info = %i") % info;
  }

  // convert the pfaffian value into mantissa/exponent form and return it
  return pcps::MantExp<double>(pfaffian[0], int(pfaffian[1])) * std::pow(double(10.0), pfaffian[1] - double(int(pfaffian[1])));

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Computes the pfaffian of a complex matrix A.
///
/// \param[in]     n        Dimension of the matrix.
/// \param[in,out] A        On entry, the n*n matrix whose pfaffian is to be computed, which is
///                         assumed to be skew-symmetric.  On exit, the matrix is overwritten.
/// \param[out]    iwork    Integer workspace, size n.
///
/// \return The pfaffian of A in mantissa/exponent format.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
pcps::MantExp<std::complex<double> > pcps::pfn::compute_pfaffian(const int n,
                                                                 std::complex<double> * const A,
                                                                 int * const iwork)
{

  // odd dimension matrices have pfaffian equal to zero
  if (n % 2 != 0)
    return pcps::MantExp<std::complex<double> >( std::complex<double>(0.0, 0.0) );

  // compute the pfaffian
  int info;
  std::complex<double> pfaffian[2];
  std::complex<double> work[2];
  double * rwork;  // not referenced for method = "P"
  int lwork = 2;
  zskpf10_("U", "P", &n, A, &n, pfaffian, iwork, work, &lwork, rwork, &info);

  // if the pfaffian fails to compute, return zero so that initial configurations using the pfaffian value will be rejected
  if (info != 0) {
    return pcps::MantExp<std::complex<double> >( std::complex<double>(0.0, 0.0) );
    //throw pcps::Exception("zskpf10_ failed in pcps::pfn::compute_pfaffian with info = %i") % info;
  }

  // check that the exponent is purely real
  if (pfaffian[1].imag() != 0.0)
    throw pcps::Exception("exponent has an imaginary component of %.2e in pcps::pfn::compute_pfaffian") % pfaffian[1].imag();

  // convert the pfaffian value into mantissa/exponent form and return it
  return pcps::MantExp<std::complex<double> >(pfaffian[0], int(pfaffian[1].real()))
          * std::pow(double(10.0), pfaffian[1] - std::complex<double>(double(int(pfaffian[1].real())), 0.0) );

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Computes the pfaffian value for a given configuration.
///
/// \param[in]     np             number of particles
/// \param[in]     ns             number of lattice sites
/// \param[in]     config         lattice configuration, size ns   
/// \param[in]     pairing_mat    the full skew-symmetric pairing matrix, size ns*ns
/// \param[out]    occ_sites      workspace for holding the occupied sites, size np   
/// \param[out]    iwork          integer workspace used in computing the pfaffian, size ns   
/// \param[out]    occ_par_mat    used to hold the occupied pairing matrix, size np*np
///
/// \return The pfaffian of the occupied pairing matrix for the given configuration
///         in mantissa/exponent format.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> pcps::MantExp<S> pcps::pfn::get_pfaffian_value(const int np,
                                                                           const int ns,
                                                                           const int * const config,
                                                                           const S * const pairing_mat,
                                                                           int * const occ_sites,            
                                                                           int * const iwork,                
                                                                           S * const occ_par_mat)       
{

  // Populate occ_sites
  pcps::pfn::init_occ_sites(np, ns, config, occ_sites, iwork);

  // build the occupied pairing matrix
  pcps::pfn::build_occ_par_mat(np, ns, pairing_mat, occ_sites, occ_par_mat);

  // compute and return the pfaffian of the occupied pairing matrix
  return pcps::pfn::compute_pfaffian(np, occ_par_mat, iwork);

}

// explicitly instantiate the function
template pcps::MantExp<double> pcps::pfn::get_pfaffian_value(const int,
                                                             const int,
                                                             const int * const,
                                                             const double * const,
                                                             int * const,
                                                             int * const,
                                                             double * const);
template pcps::MantExp<std::complex<double> > pcps::pfn::get_pfaffian_value(const int,
                                                                            const int,
                                                                            const int * const,
                                                                            const std::complex<double> * const,
                                                                            int * const,
                                                                            int * const,
                                                                            std::complex<double> * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Computes the ratio of the pfaffians for two lattice configurations that differ
///          at the specified sites.
///
/// \param[in]     np             number of particles
/// \param[in]     ns             number of lattice sites
/// \param[in]     nd             the number of sites that differ between the configurations
/// \param[in]     diff_sites     the indices of the differing sites ( 0 <= index < ns ), size nd
/// \param[in]     new_config     new lattice configuration, size ns
/// \param[in]     old_config     old lattice configuration, size ns
/// \param[in]     pairing_mat    the full skew-symmetric pairing matrix, size ns*ns
/// \param[in,out] inv_mat        The inverse of the occupied pairing matrix for old_config,
///                               size np*np.  
///                               Upon exit, the rows/columns of this matrix may be reordered.
/// \param[in,out] occ_sites      The list of occupied sites, size np.
///                               Upon exit, the order of sites in this list may change to keep
///                               track of row/column reorderings in inv_mat.
/// \param[in,out] row_indices    The row indices for occupied sites, size ns.
///                               Upon exit, this list is updated to reflect
///                               any row/column reorderings in inv_mat.
/// \param[out]    wv1            Work space, size np.
///
/// \return The ratio of the pfaffians of the occupied pairing matrices at the new and old
///         configurations ( i.e. pfn(new) / pfn(old) ).
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S pcps::pfn::pfn_ratio(const int np,
                                                    const int ns,
                                                    const int nd,
                                                    const int * const diff_sites,
                                                    const int * const new_config,
                                                    const int * const old_config,
                                                    const S * const pairing_mat,
                                                    S * const inv_mat,
                                                    int * const occ_sites,
                                                    int * const row_indices,
                                                    S * const wv1)
{

  // check that the number of sites affected by the move is supported
  assert(nd == 0 || nd == 2 );

  // return one if no sites were changed
  if (nd == 0)
    return pcps::unity<S>();

  // determine which among the affected lattice sites are newly occupied, and which were previously occupied
  int new_occ_sites[2];
  int old_occ_sites[2];
  pcps::pfn::get_diff_occs(np, ns, nd, diff_sites, new_config, old_config, new_occ_sites, old_occ_sites);

  // exchange rows and columns so that the altered row and column are the first row and column
  pcps::pfn::make_first_row_col(np, ns, row_indices[old_occ_sites[0]], inv_mat, occ_sites, row_indices, wv1);

  // compute the pfaffian ratio
  S retval = pcps::unity<S>();
  for (int i = 1; i < np; i++)
    retval += inv_mat[i*np+0] * ( pairing_mat[new_occ_sites[0] * ns + occ_sites[i]] - pairing_mat[occ_sites[0] * ns + occ_sites[i]] );

  // return the pfaffian ratio
  return retval;

}

// explicitly instantiate the function
template double pcps::pfn::pfn_ratio(const int np,
                                     const int ns,
                                     const int nd,
                                     const int * const diff_sites,
                                     const int * const new_config,
                                     const int * const old_config,
                                     const double * const pairing_mat,
                                     double * const inv_mat,
                                     int * const occ_sites,
                                     int * const row_indices,
                                     double * const wv1);
template std::complex<double> pcps::pfn::pfn_ratio(const int np,
                                                   const int ns,
                                                   const int nd,
                                                   const int * const diff_sites,
                                                   const int * const new_config,
                                                   const int * const old_config,
                                                   const std::complex<double> * const pairing_mat,
                                                   std::complex<double> * const inv_mat,
                                                   int * const occ_sites,
                                                   int * const row_indices,
                                                   std::complex<double> * const wv1);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Updates the inverse of the occupied pairing matrix and the lists of which sites are
///          occupied and which occupied sites correspond to which rows upon changing to a new
///          lattice configuration.
///
/// \param[in]     np           number of particles
/// \param[in]     ns           number of lattice sites
/// \param[in]     nd           the number of sites that differ between the configurations
/// \param[in]     diff_sites   the indices of the differing sites ( 0 <= index < ns ), size nd
/// \param[in]     new_config   new lattice configuration, size ns
/// \param[in]     old_config   old lattice configuration, size ns
/// \param[in]     pairing_mat  the full skew-symmetric pairing matrix, size ns*ns
/// \param[in,out] inv_mat      Size np*np.
///                             On entry, the inverse of the occupied pairing matrix for old_config.
///                             On exit, the inverse of the occupied pairing matrix for new_config.
/// \param[in,out] occ_sites    Size np.
///                             On entry, the list of occupied sites for old_config.
///                             On exit, the list of occupied sites for new_config.
/// \param[in,out] row_indices  Size ns.
///                             On entry, the occupied sites' row indices for old_config.
///                             On exit, the occupied sites' row indices for new_config.
/// \param[out]    wv1          Work space, size np.
/// \param[out]    wv2          Work space, size np.
/// \param[out]    wv3          Work space, size np.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::pfn::update(const int np,
                                               const int ns,
                                               const int nd,
                                               const int * const diff_sites,
                                               const int * const new_config,
                                               const int * const old_config,
                                               const S * const pairing_mat,
                                               S * const inv_mat,
                                               int * const occ_sites,
                                               int * const row_indices,
                                               S * const wv1,
                                               S * const wv2,
                                               S * const wv3)
{

  // check that the number of sites affected by the move is supported
  assert(nd == 0 || nd == 2 );

  // nothing to do if no sites are affected
  if (nd == 0) return;

  // determine which among the affected lattice sites are newly occupied, and which were previously occupied
  int new_occ_sites[2];
  int old_occ_sites[2];
  pcps::pfn::get_diff_occs(np, ns, nd, diff_sites, new_config, old_config, new_occ_sites, old_occ_sites);

  // exchange rows and columns so that the altered row and column are the first row and column
  pcps::pfn::make_first_row_col(np, ns, row_indices[old_occ_sites[0]], inv_mat, occ_sites, row_indices, wv1);

  // get row difference (which is the negative of the column difference due to antisymmetry of the occupied pairing matrix)
  wv1[0] = pcps::zero<S>();
  for (int i = 1; i < np; i++)
    wv1[i] = pairing_mat[new_occ_sites[0] * ns + occ_sites[i]] - pairing_mat[occ_sites[0] * ns + occ_sites[i]];

  // product of inverse with row difference
  pcps::xgemm('N', 'N', np, 1, np, pcps::unity<S>(), inv_mat, np, wv1, np, pcps::zero<S>(), wv2, np);

  // get first column of the inverse matrix
  pcps::xcopy(np, inv_mat, np, wv3, 1);

  // update after row change
  pcps::xgemm('N', 'N', np, np, 1, -1.0 / ( pcps::unity<S>() + wv2[0] ), wv2, np, wv3, 1, pcps::unity<S>(), inv_mat, np);

  // product of inverse with column difference
  pcps::xgemm('T', 'N', np, 1, np, -pcps::unity<S>(), inv_mat, np, wv1, np, pcps::zero<S>(), wv2, np);

  // get first row of the inverse matrix
  pcps::xcopy(np, inv_mat, 1, wv3, 1);

  // update after row change
  pcps::xgemm('N', 'N', np, np, 1, -1.0 / ( pcps::unity<S>() + wv2[0] ), wv3, np, wv2, 1, pcps::unity<S>(), inv_mat, np);

  // update occ_sites and row_indices
  row_indices[old_occ_sites[0]] = -1;
  row_indices[new_occ_sites[0]] =  0;
  occ_sites[0] = new_occ_sites[0];

}

// explicitly instantiate the function
template void pcps::pfn::update(const int np,
                                const int ns,
                                const int nd,
                                const int * const diff_sites,
                                const int * const new_config,
                                const int * const old_config,
                                const double * const pairing_mat,
                                double * const inv_mat,
                                int * const occ_sites,
                                int * const row_indices,
                                double * const wv1,
                                double * const wv2,
                                double * const wv3);
template void pcps::pfn::update(const int np,
                                const int ns,
                                const int nd,
                                const int * const diff_sites,
                                const int * const new_config,
                                const int * const old_config,
                                const std::complex<double> * const pairing_mat,
                                std::complex<double> * const inv_mat,
                                int * const occ_sites,
                                int * const row_indices,
                                std::complex<double> * const wv1,
                                std::complex<double> * const wv2,
                                std::complex<double> * const wv3);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Computes the ratios of the pfaffian's derivatives with respect to the pairing
///          matrix elements to the pfaffian's value.
///
/// \param[in]     np           number of particles
/// \param[in]     ns           number of lattice sites
/// \param[in]     inv_mat      The inverse of the occupied pairing matrix for the current
///                             lattice configuration, size np*np.
/// \param[in]     occ_sites    The list of occupied sites for the current lattice configuration,
///                             size np.
/// \param[in]     row_indices  The occupied sites' row indices for the current lattice
///                             configuration, size ns.
/// \param[out]    ratios       On exit, the derivative ratios, size ns*(ns-1)/2.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::pfn::der_ratios(const int np,
                                                   const int ns,
                                                   const S * const inv_mat,
                                                   const int * const occ_sites,
                                                   const int * const row_indices,
                                                   S * const ratios)
{

  // zero the derivatives
  pcps::xscal(ns*(ns-1)/2, pcps::zero<S>(), ratios, 1);

  // add the derivative ratios of the occupied pairing matrix's pfaffian
  for (int i = 0; i < np; i++)
  for (int j = 0; j < np; j++) {
    const int I = occ_sites[i];
    const int J = occ_sites[j];
    if (I < J)
      ratios[pcps::cmpd_pair_index(I,J,ns)] += 0.5 * ( inv_mat[j*np+i] - inv_mat[i*np+j] );
  }

}

// explicitly instantiate the function
template void pcps::pfn::der_ratios(const int np,
                                    const int ns,
                                    const double * const inv_mat,
                                    const int * const occ_sites,
                                    const int * const row_indices,
                                    double * const ratios);
template void pcps::pfn::der_ratios(const int np,
                                    const int ns,
                                    const std::complex<double> * const inv_mat,
                                    const int * const occ_sites,
                                    const int * const row_indices,
                                    std::complex<double> * const ratios);
